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Abstract — We  propose  algorithms  for  nonparametric  sample- 
based  spacial  change  detection  and  estimation  in  large  scale 
sensor  networks.  We  collect  random  samples  containing  the 
location  of  sensors  and  their  local  decisions,  and  assume  that 
the  local  decisions  can  be  “stimulated”  or  “normal”,  reflecting 
the  local  strength  of  some  stimulating  agent.  Then  change  in 
the  location  of  the  agent  manifests  itself  by  a  change  in  the 
distribution  of  “stimulated”  sensors.  In  this  paper,  we  are  aiming 
at  developing  a  test  that,  given  two  collections  of  samples, 
can  decide  whether  the  distribution  generating  the  samples  has 
changed  or  not,  and  give  an  estimated  changed  area  if  a  change 
is  indeed  detected.  The  focus  of  this  paper  is  to  reduce  the 
complexity  of  the  detection  and  estimation  algorithm.  We  propose 
two  fast  algorithms  with  almost  linear  complexity  and  analyze 
their  completeness,  flexibility  and  robustness. 

Keywords:  Non-parametric  statistical  methods,  Change  detec¬ 
tion,  Spacial  change,  Sensor  Networks. 


I.  Introduction 

Given  a  large  scale  sensor  network,  where  sensors  are 
distributed  to  detect  the  local  presence  of  some  agents,  we 
are  interested  in  deciding  whether  the  location  of  agents  has 
changed  by  looking  at  random  samples  drawn  from  sensors. 
We  assume  that  the  agents  have  local  property,  i.e.  their  most 
intensive  areas  are  clustered  in  space,  forming  some  “heated 
areas”  in  the  sensor  field.  One  example  is  the  presence  of 
targets.  Another  example  is  sources  stimulating  “excited  state” 
of  sensors,  with  this  stimulating  effect  decaying  as  distance 
grows.  A  sensor  reports  “RED”  if  it  detects  the  agent  or 
“GREEN”  otherwise.  Using  the  architecture  of  SENMA  [6], 
and  either  letting  sensors  transmitting  their  location  infor¬ 
mation  together  with  decisions  or  polling  sensors  at  specific 
locations,  we  can  collect  a  set  of  samples  containing  local 
decisions  and  the  locations  where  these  decisions  are  made. 

The  problem  of  change  detection  and  estimation  in  SENMA, 
illustrated  in  Figure  1,  is  to  decide,  from  two  consecutive 
records,  whether  there  is  a  change  of  the  underlying  probabil¬ 
ity  distribution  (detection)  and  the  distribution  of  the  change 
(estimation)  if  changes  occur.  In  practice,  random  access 

This  work  is  supported  by  the  U.  S.  Army  Research  Laboratory  under 
the  Collaborative  Technology  Alliance  Program,  Cooperative  Agreement 
DAAD 19-01  -2-001 1.  The  U.  S.  Government  is  authorized  to  reproduce  and 
distribute  reprints  for  Government  purposes  notwithstanding  any  copyright 
notation  thereon. 


First  data  collection  Second  data  collection 

Fig.  1.  Sensor  Networks  with  Mobile  Access.  RED:  sensors  with  detection. 
GREEN:  sensors  with  no  detection.  Gray:  Data  not  collected. 

is  used  at  the  mobile  access  point  when  collecting  packets 
from  the  sensor  field(See  Fig.  1),  and  sensors  can  make 
mistakes  in  their  decisions,  i.e.  reporting  “RED”  with  no  agent 
present  (false  alarm)  or  “GREEN”  when  there  is  an  agent 
present  (missing  detection).  Hence  neither  is  there  guarantee 
that  a  measurement  at  a  particular  sensor  in  one  data  collection 
will  appear  in  the  next,  nor  can  we  say  a  certain  measurement 
is  true  or  not.  The  remedy  is  to  take  a  large  number  of 
samples  to  compensate  for  the  ill  effect  of  different  locations 
and  false  reports.  However,  two  factors  must  be  considered 
for  change  detection  and  estimation  in  a  large  scale  sensor 
networks:  the  time  required  or  the  number  of  packets  required 
for  data  collection  and  the  complexity  of  the  detection  and 
estimation  algorithm.  In  practice,  minimizing  the  number  of 
samples  required  for  detection  reduces  the  collection  time 
of  a  mobile  AP  and,  more  importantly  in  terms  of  energy 
consumption,  the  number  of  transmissions  from  sensors.  In 
a  previous  paper  [1],  we  have  provided  a  mathematical  char¬ 
acterization  of  the  required  sample  size  for  which,  with  high 
probability,  the  distance  of  any  two  probability  distributions 
can  be  estimated  accurately  from  empirical  distance  between 
samples  generated  by  these  distributions.  The  basic  tool  used 
here  is  the  Vapnik-Chervonenkis  Theory  [2],  In  this  paper,  we 
are  mainly  concerned  with  the  complexity  of  the  detection  and 
estimation  algorithm.  We  present  variations  of  an  algorithm 
in  [1]  and  analyze  their  complexity,  flexibility  and  levels  of 
confidence. 

For  many  applications,  there  is  no  prior  knowledge  about 
the  distribution  of  the  sensor  states.  The  change  detection 
and  estimation  is  therefore  non-parametric  in  which  we  make 
no  assumption  about  the  specific  form  of  the  probability 
distribution  of  the  binary  random  field.  See  the  introduction  of 
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[1]  for  a  brief  review  of  some  previous  work  in  non-parametric 
change  detection  problem. 


II.  Algorithms 


In  this  section,  we  take  one  of  the  three  algorithms  in  [1] 
and  derive  its  two  variations  with  marginal  performance  but 
greatly-reduced  complexity.  The  algorithm  we  are  based  on  is 
Algorithm  3:  Search  in  Axis-aligned  Rectangles.  Using  the 
same  boundary  search  idea  in  Algorithm  3,  we  reduce  its 
complexity  to  be  almost  linear  by  projecting  samples  onto 
one-dimension,  and  refining  the  results  by  multiple  projections 
at  a  cost  of  only  a  constant  factor  increase  in  running  time. 
As  algorithms  in  [1],  both  of  the  algorithms  are  based  on  the 
fact  that  empirical  probability  distributions  of  i.i.d.  samples 
uniformly  converge  to  their  actual  distributions,  no  matter  what 
kind  of  distributions  they  are  [2],  Specifically,  we  rely  on 
results  in  [1]  to  provide  a  guaranteed  level  of  accuracy. 

Practical  implementations  require  that  the  distance  between 
two  empirical  distributions  be  computed  in  a  timely  manner. 
Given  limited  computation  resource,  it  is  of  interest  to  look 
at  the  largest  sample  size  we  can  handle  with  our  detec¬ 
tion/estimation  algorithms.  We  therefore  need  to  reduce  the 
complexity  of  our  algorithms  in  terms  of  sample  size  to  enable 
efficient  processing  of  large  samples.  The  key  to  complexity 
reduction  is  the  projection  of  samples  onto  one-dimension. 
After  the  projection,  the  search  is  reduced  to  search  on  a  line, 
and  by  reusing  previous  computation  the  search  is  done  in  a 
linear  number  of  steps. 

We  consider  the  following  detector: 

Given  two  collection  of  samples  Si  and  S2,  drawn  i.i.d  from 
probability  distributions  P\  and  P-2  respectively,  and  threshold 
e  €  (0, 1),  for  simple  binary  hypothesis  Ho  :  P\  =  P2  vs. 
Hi  :  Pi  /  P2,  the  detector  is  defined  as 

<5 (Si,  S'2;  e)  =  1  if  diffmax  >  e 

=  0  otherwise 

where  diffmax  is  the  maximum  difference  in  the  empirical 
probabilities  computed  by  some  algorithm. 

Furthermore,  if  5 (Si,  S'2;  e)  =  1,  then  the  estimated  changed 
area  is  the  area  where  diffmax  is  achieved. 

Now  the  problem  is  reduced  to  finding  diffmax.  Previous 
algorithmsfsee  [1])  include  Algorithm  1:  Exhaustive  Search 
in  Planar  Disks.  Algorithm  2:  Search  in  Sample-centered 
Disks  and  Algorithm  3:  Search  in  Axis-aligned  Rectangles, 
with  complexities  0(M4 ),  0(M2  log  M)  and  0(M3),  where 
M  =  |  Si  U  S2I.  Based  on  the  same  boundary  search  idea  of 
Algorithm  3,  we  have  the  following  faster  algorithms. 


Algorithm  3.1:  Search  in  Axis-aligned  Stripes 

This  is  a  simplified  version  of  Algorithm  3.  The  basic  idea 
is  to  project  samples  onto  x  and  y  coordinates,  and  perform 
change  detection/estimation  on  each  coordinate.  The  search  is 
based  on  the  same  boundary  search  idea  of  Algorithm  3,  but 
by  projection  the  complexity  of  the  search  is  reduced  to  be 
linear. 

Let  A  be  the  collection  of  vertical  stripes,  i.e.,  axis  aligned 
rectangles  with  height  equal  to  the  field  height.  Similarly,  let 
B  be  the  collection  of  horizonal  stripes.  Given  a  collection 
of  samples  S  =  Si  U  S'2,  it  suffices  to  consider  finite  subsets 


H%(S)  C  A  and  H\(S)  C  B  defined  by 

HUS) 

{ V (x*i ,  Xj )  .  si  (xi,yf), 

si  =  (xPVj)  e  S} 

(1) 

HUS) 

=  {H(yk,yi)  ■■  Sk  =  (xk,yk), 

Si  =  ( Xi,yi )  e  S} 

(2) 

where  V(xi,xf) 

is  the  vertical  stripe  with  left 

and  right 

boundary  Xi  and  Xj,  and  H(yk,  y{)  is  the  horizonal  stripe  with 
lower  and  upper  boundary  yk  and  ly.  See  Figure  2. 


Fig.  2.  Members  of  and  7-f| 

Recall  the  definition  of  completeness  in  [1],  we  have  the 
following  claim. 

Claim  0.1:  Let  A.  B  be  as  defined  above.  Given  Si,  S'2,  the 
finite  subset  H^(Si  U  S2)  U  7f|(Si  U  S2)  defined  in  1  and  2 
is  complete  for  A  U  B  w.r.t.  Si  U  S2. 

Given  S  =  Si  U  S2,  Algorithm  3.1  performs  the  following 
exhaustive  search  in  H%(S)  U  H\(S)\ 

Si  (T  A\  |S2nA| 

~W\  \sT  ' 

where  Si,  S2  are  sets  of  the  previous  and  current  red  samples. 

The  algorithm  includes  (i)  projecting  the  red  samples  onto 
x  and  y  coordinates;  (ii)  sorting  the  projected  samples  into 
increasing  order;  (iii)  in  x  coordinate,  scanning  from  left  to 
right,  at  the  ith  sample  with  x  coordinate  x^  computing  Fx(i), 


max 


the  difference  in  the  empirical  probability  of  stripe  V(0,  &*) 
between  Si  and  S2,  which  is  given  by 

Fx(i)  =  Fx(i  -  1)  +  T^-r  if  Sies1  (3) 
l*i  I 

=  Fx(i  -  1)  -  T^-r  if  Si  e  S2  (4) 

1*2 1 

computing  Fy(-)  in  y  coordinate  similarly;  (iv)  finding  the 
peak  and  bottom  point  of  Fx(-)  achieved  at  mi  and  m2,  and 
the  peak  and  bottom  of  Fy(-)  at  m  and  77.2.  Then  the  max¬ 
imum  difference  diffmax  =  max^gH«uw* 
is  given  by 

max(Fx(mi)  -  Fx(m2),  Fy(n  1)  -  Fy(n2 )) 

and  the  axis  aligned  stripe  bounded  by  x  =  xmi  and  x  =  x„,2 
if  Fx(mi)  -  Fx(m2 )  >  ^(m)  -  Fy(n2),  or  y  =  yni  and 
V  =  yn 2  if  Fx(mi)  -  Fx(m2)  <  Fy(m)  -  Fv(n2),  gives  an 
estimation  of  the  changed  area. 

The  complexity  of  Algorithm  3.1  for  sample  size  M  = 
| S1!  US2|  is  O(MlogM)1.  Thus  by  projection  we  reduce  the 
complexity  to  make  Algorithm  3.1  an  almost-linear  algorithm. 
Like  Algorithm  3,  Algorithm  3.1  is  also  not  amendable  to 
relativized  discrepancy  defined  in  [1]. 

Algorithm  3.2;  Search  in  Random  Stripes 

This  is  a  variation  of  Algorithm  3.1.  It  is  based  on  the 
same  projection  and  boundary  search  idea  as  in  Algorithm 
3.1.  The  difference  is  when  performing  the  projection,  we 
project  samples  onto  a  random  direction  instead  of  the  fixed 
direction  of  x  or  y  axis.  The  rest  of  the  algorithm  is  the  same 
as  Algorithm  3.1. 

Algorithm  3.2  has  the  same  complexity  and  flexibility  as 
Algorithm  3.1.  It  is  also  complete  for  the  class  of  stripes 
along  the  particular  direction  as  chosen  by  Algorithm  3.2 
in  the  sense  that  it  exhausts  all  such  stripes  with  at  least 
one  sample  point  on  each  boundary.  Its  advantage  is  that  it 
is  robust  with  respect  to  changing  patterns.  This  is  in  the 
sense  that  Algorithm  3.2  will  perform  equally  well  under  all 
kinds  of  changing  patterns(the  directions  along  which  changes 
occur)  while  Algorithm  3.1  can  be  affected  significantly  by 
the  particular  changing  pattern  that  is  taken.  For  example. 
Algorithm  3.1  is  vulnerable  to  the  pattern  where  changes 
always  occur  along  a  tilted  line  of  45°  or  135°,  because  in 
that  case  the  increasing  and  decreasing  parts  of  the  change 
will  largely  get  cancelled  when  projected  onto  axes.  We  can 
also  project  onto  multiple  random  directions  to  increase  the 
accuracy  of  the  algorithm  at  the  cost  of  a  constant  factor 
increase  in  the  complexity.  For  example,  in  the  simulation  we 

'The  complexity  0(M  log  M)  is  determined  by  the  critical  step  of  sorting. 
Strictly  speaking,  this  is  more  than  linear  complexity.  However,  for  sample  size 
of  practical  interest,  the  factor  log  M  can  be  incorporated  into  the  constant 
factor  and  the  algorithm  scales  like  a  linear  algorithm. 


use  two  random  directions  for  fair  comparison  with  Algorithm 
3.1  since  the  latter  use  two  directions,  i.e.  x  and  y  axes. 


III.  Simulation 

A.  Simulation  Setup 

In  the  simulation,  we  assume  that  the  sensor  field  is  a  50  x  50 
to2  square,  and  sensors  are  randomly  distributed  in  the  field. 
Assume  the  distribution  of  operating  sensors  does  not  change. 
Then  it  suffices  to  consider  only  “RED”  samples.  We  adopt 
a  simple  sampling  strategy  that  is  continuing  sampling  until 
enough  “RED”  samples  are  collected.  The  change  simulated 
here  is  a  change  in  the  location  of  the  most  intensive  area 
of  the  stimulating  agent.  For  the  stimulating  agent,  we  adopt 
a  simplified  model  that  says  a  sensor  reports  “RED”  with 
probability  p  if  it  falls  inside  distance  r  from  the  center  of 
the  agent  and  q  if  it  falls  outside  this  distance,  where  p  >  q. 
We  let  p  <  1,  q  >  0  to  model  the  randomness  in  the  decision 
of  each  sensor. 

Given  the  highest  false  alarm  allowed  (size  of  the  detector), 
we  first  perform  Monte  Carlo  runs  to  decide  the  detection 
threshold,  and  then  do  detection  and  estimation  using  this 
threshold.  Hence  the  nonparametric  property  of  the  simulation 
is  preserved.  Difference  between  algorithms  becomes  in  which 
class(e.g.  disks,  stripes,  rectangles)  we  are  trying  to  fit  the 
changed  area  and  how  to  find  the  best  match  in  the  class. 

In  our  simulation,  we  let  p  =  0.8,  q  =  0.2,  r  =  ^  m.  The 
size  of  the  detector  is  a  =  0.05. 


B.  Detection  Threshold  and  Missing  Detection  Probability 

In  the  simulation,  we  first  fix  a  range  of  sample  size  we 
are  interested  in,  and  then  do  Monte  Carlo  runs  at  each 
sample  size  with  fixed  agent  location  (no  change)  to  find  out 
the  distribution  of  dijfmaxi the  maximum  change  in  empirical 
probability  found  by  our  algorithms)  under  'Ho-  The  detection 
threshold  for  that  sample  size  is  determined  as  the  smallest 
value  such  that  the  frequency  that  diffmax  exceeds  this  value 
is  below  the  detector’s  size  a. 

For  evaluation  purpose,  we  compare  the  detection  thresholds 
and  missing  detection  probabilities  for  Algorithm  3.1  and 
Algorithm  3.2.  Specifically,  we  plot  the  detection  threshold 
obtained  by  simulation  as  described  above  as  a  function  of 
the  red  sample  size  for  red  sample  size  between  594  and  104, 
and  compare  it  to  the  theoretical  threshold  computed  from  the 
upper  bound  in  [1],  In  computing  the  theoretical  threshold,  we 
find  the  smallest  e  such  that  8(2rc)2e_ne  Z32  <  a  and  choose 
e  as  the  theoretical  threshold.  See  Fig.  3. 
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Detection  Threshold  vs.  Sample  Size  false  alarm:  a  =  0.05 
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Fig.  3.  Detection  threshold  as  a  function  of  the  red  sample  size 


As  expected,  both  simulated  and  theoretical  threshold  is  a 
decreasing  function  of  the  sample  size,  which  reflects  a  trade¬ 
off  between  sampling  time,  energy  consumed  and  data  pro¬ 
cessing  expense  and  detection  precision.  Also  as  expected,  the 
theoretical  threshold  is  a  loose  upper  bound  of  the  simulated 
threshold.  This  is  because  of  the  nonparametric  nature  of  the 
theoretical  upper  bound.  This  bound  is  claimed  to  hold  under 
arbitrary  distributions  by  the  Vapnik-Chervonenkis  Theory. 
Therefore  for  a  given  distribution,  this  bound  is  loose. 


It  is  worth  noting  that  here  Algorithm  3.1  performs  slightly 
better  because  of  the  particular  set  up  we  are  using  in 
the  simulation.  For  randomly  changing  locations.  Algorithm 
3.1  ensures  good  utilization  of  samples  by  searching  along 
orthogonal  directions,  while  Algorithm  3.2  performs  equally 
well  under  all  changing  directions  by  searching  along  random 
directions. 
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Then  using  the  threshold  obtained  from  simulation,  for 
randomly  changing  locations,  we  see  from  Fig.  4  how  the 
missing  detection  probability  of  these  algorithms  decreases  as 
sample  size  grows. 


Missing  Detection  vs.  Sample  Size  false  alarm:  a  =  0.05 
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Fig.  4.  Missing  detection  probability  as  a  function  of  the  red  sample  size 


